pacman::p_load(sf, raster, spatstat, tmap, tidyverse, dplyr, readr, spNetwork)Take-Home Ex01
Geospatial Analytics for Social Good: Thailand Road Accident Case Study
1.0 Installing and Loading the R packages
2.0 Importing the spatial data
2.1 Importing in Thailand road accidents in Bangkok province between 2019-2022
The code below filters out for NA and blanks in the longitude, latitude, and province columns. It also filters specifically for “Bangkok” province, which is the area of this study. The coordinate reference system is specified as EPSG:32647 (UTM zone 47N coordinate system), which is suitable for spatial analyses within the Bangkok region. This is because WGS84 is not ideal for precise distance and area calculations as it represents Earth as a curved surface, UTM, on the other hand, projects the Earth onto a flat plane, which reduces distortions and allows for more accurate measurements. also, with mutate function, we can get month and day of week of the Bangkok accident #label = true make month factor, abbr= abbrevation
rdacc_sf <- read_csv("C:/jhui-chen/ISSS626-GAA/Take-Home_Ex/Take-Home_Ex01/data/thai_road_accident_2019_2022.csv") %>%
mutate(month = month(incident_datetime)) %>%
mutate(monthfac = month(incident_datetime, label = TRUE, abbr = TRUE)) %>%
mutate(dayofweek = day(incident_datetime)) %>%
filter(!is.na(longitude) & longitude != "", # filter out NA and blanks in longitude
!is.na(latitude) & latitude != "", # filter out NA and blanks in latitude
!is.na(province_en) & province_en != "", # filter out NA and blanks in province
province_en == "Bangkok") %>% # filter for province "Bangkok"
st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>% #convert data frame into simple feature (sf). Tells R which columns are long and lat, and crs = 4326 specifies the coordinate reference system are in WGS 84.
st_transform(crs = 32647) #projects spatial data into the UTM zone 47N coordinate system.Rows: 81735 Columns: 18
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (10): province_th, province_en, agency, route, vehicle_type, presumed_c...
dbl (6): acc_code, number_of_vehicles_involved, number_of_fatalities, numb...
dttm (2): incident_datetime, report_datetime
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#view(rdacc_sf)Checks on coordinate reference system.
st_crs(rdacc_sf) Coordinate Reference System:
User input: EPSG:32647
wkt:
PROJCRS["WGS 84 / UTM zone 47N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 47N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",99,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand."],
BBOX[0,96,84,102]],
ID["EPSG",32647]]Having a look.
rdacc_sfSimple feature collection with 6089 features and 19 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 644139 ymin: 1505014 xmax: 708461.2 ymax: 1542968
Projected CRS: WGS 84 / UTM zone 47N
# A tibble: 6,089 × 20
acc_code incident_datetime report_datetime province_th province_en
* <dbl> <dttm> <dttm> <chr> <chr>
1 629691 2019-01-01 03:05:00 2019-01-01 03:05:00 กรุงเทพมหานคร Bangkok
2 629689 2019-01-01 05:42:00 2019-01-01 05:42:00 กรุงเทพมหานคร Bangkok
3 604307 2019-01-01 10:10:00 2019-03-06 10:49:00 กรุงเทพมหานคร Bangkok
4 3793736 2019-01-01 17:30:00 2020-03-11 13:15:00 กรุงเทพมหานคร Bangkok
5 599070 2019-01-01 19:20:00 2019-01-01 20:57:00 กรุงเทพมหานคร Bangkok
6 613605 2019-01-01 21:40:00 2019-11-18 10:57:00 กรุงเทพมหานคร Bangkok
7 615219 2019-01-01 22:15:00 2019-12-27 10:12:00 กรุงเทพมหานคร Bangkok
8 613557 2019-01-01 23:05:00 2019-11-18 10:57:00 กรุงเทพมหานคร Bangkok
9 629699 2019-01-02 03:40:00 2019-01-02 03:40:00 กรุงเทพมหานคร Bangkok
10 3793707 2019-01-02 07:00:00 2020-03-11 12:57:00 กรุงเทพมหานคร Bangkok
# ℹ 6,079 more rows
# ℹ 15 more variables: agency <chr>, route <chr>, vehicle_type <chr>,
# presumed_cause <chr>, accident_type <chr>,
# number_of_vehicles_involved <dbl>, number_of_fatalities <dbl>,
# number_of_injuries <dbl>, weather_condition <chr>, road_description <chr>,
# slope_description <chr>, month <dbl>, monthfac <ord>, dayofweek <int>,
# geometry <POINT [m]>plot(st_geometry(rdacc_sf))
2.2 Importing in Thailand road shapefile
Importing the Thailand road shapefile The coordinate reference system is specified as EPSG:32647 (UTM zone 47N coordinate system), which is suitable for spatial analyses within the Bangkok region.
road <- st_read(dsn = "C:/jhui-chen/ISSS626-GAA/Take-Home_Ex/Take-Home_Ex01/data",
layer = "hotosm_tha_roads_lines_shp")Reading layer `hotosm_tha_roads_lines_shp' from data source
`C:\jhui-chen\ISSS626-GAA\Take-Home_Ex\Take-Home_Ex01\data'
using driver `ESRI Shapefile'
Simple feature collection with 2792590 features and 14 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: 97.34457 ymin: 5.643645 xmax: 105.6528 ymax: 20.47168
CRS: NAst_crs(road) <- 4326 #Set the CRS manually to WGS 84 (EPSG:4326)
road <- st_transform(road, crs = 32647) #Transform the CRS to UTM zone 47N (EPSG:32647)Checks on coordinate reference system.
st_crs(road)Coordinate Reference System:
User input: EPSG:32647
wkt:
PROJCRS["WGS 84 / UTM zone 47N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 47N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",99,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand."],
BBOX[0,96,84,102]],
ID["EPSG",32647]]Having a look.
roadSimple feature collection with 2792590 features and 14 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: 325313.7 ymin: 624248.4 xmax: 1215576 ymax: 2263968
Projected CRS: WGS 84 / UTM zone 47N
First 10 features:
name name_en highway surface smoothness
1 ถนนฉลองกรุง Chalong Krung Road secondary paved <NA>
2 ซอยฉลองกรุง 1/1 Soi Chalong Krung 1/1 residential <NA> <NA>
3 <NA> <NA> secondary_link <NA> <NA>
4 <NA> <NA> service <NA> <NA>
5 ถนนฉลองกรุง Chalong Krung Road secondary concrete <NA>
6 <NA> <NA> service <NA> <NA>
7 ถนนเอราวัณ 1 Erawan 1 Road tertiary <NA> <NA>
8 <NA> <NA> path unpaved <NA>
9 <NA> <NA> service <NA> <NA>
10 <NA> <NA> residential <NA> <NA>
width lanes oneway bridge layer source name_th osm_id osm_type
1 <NA> <NA> yes <NA> <NA> <NA> ถนนฉลองกรุง 1125681229 ways_line
2 <NA> <NA> <NA> <NA> <NA> <NA> ซอยฉลองกรุง 1/1 594401607 ways_line
3 <NA> <NA> yes <NA> <NA> <NA> <NA> 472283206 ways_line
4 <NA> <NA> <NA> <NA> <NA> <NA> <NA> 594401608 ways_line
5 <NA> 2 yes yes 1 Bing ถนนฉลองกรุง 116847248 ways_line
6 <NA> <NA> <NA> <NA> <NA> <NA> <NA> 317485095 ways_line
7 <NA> <NA> <NA> <NA> <NA> <NA> ถนนเอราวัณ 1 378672881 ways_line
8 <NA> <NA> <NA> <NA> <NA> GPS <NA> 1238351123 ways_line
9 <NA> <NA> <NA> <NA> <NA> <NA> <NA> 909942692 ways_line
10 <NA> <NA> <NA> <NA> <NA> <NA> <NA> 694824299 ways_line
geometry
1 MULTILINESTRING ((693686.1 ...
2 MULTILINESTRING ((693358 15...
3 MULTILINESTRING ((692949.1 ...
4 MULTILINESTRING ((693256 15...
5 MULTILINESTRING ((692810.8 ...
6 MULTILINESTRING ((693877.2 ...
7 MULTILINESTRING ((677182.3 ...
8 MULTILINESTRING ((486572.6 ...
9 MULTILINESTRING ((629009.2 ...
10 MULTILINESTRING ((629703.9 ...#plot(road)2.3 Importing in Thailand - Subnational Administrative Boundaries
Investigating the layers available in my local drive. The “tha_admbnda_adm1_rtsd_20220121” layer has 77 features, most likely being the provincial boundaries, as there are 77 provinces in Thailand and Bangkok is one of them.
gdb_layers <- st_layers(dsn = "C:/jhui-chen/ISSS626-GAA/Take-Home_Ex/Take-Home_Ex01/data")
print(gdb_layers)Driver: ESRI Shapefile
Available layers:
layer_name geometry_type features fields crs_name
1 hotosm_tha_roads_lines_shp Line String 2792590 14 <NA>
2 tha_admbnda_adm0_rtsd_20220121 Polygon 1 13 WGS 84
3 tha_admbnda_adm1_rtsd_20220121 Polygon 77 16 WGS 84
4 tha_admbnda_adm2_rtsd_20220121 Polygon 928 19 WGS 84
5 tha_admbnda_adm3_rtsd_20220121 Polygon 7425 22 WGS 84
6 tha_admbndl_admALL_rtsd_itos_20220121 Line String 22425 5 WGS 84
7 tha_admbndp_admALL_rtsd_itos_20220121 Point 7425 23 WGS 84
8 tha_admbndt_adminUnitLookup NA 4 6 <NA>province <- st_read(dsn = "C:/jhui-chen/ISSS626-GAA/Take-Home_Ex/Take-Home_Ex01/data", layer = "tha_admbnda_adm1_rtsd_20220121")Reading layer `tha_admbnda_adm1_rtsd_20220121' from data source
`C:\jhui-chen\ISSS626-GAA\Take-Home_Ex\Take-Home_Ex01\data'
using driver `ESRI Shapefile'
Simple feature collection with 77 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 97.34336 ymin: 5.613038 xmax: 105.637 ymax: 20.46507
Geodetic CRS: WGS 84st_crs(province) <- 4326 #Set the CRS manually to WGS 84 (EPSG:4326)
province <- st_transform(province, crs = 32647) #Transform the CRS to UTM zone 47N (EPSG:32647)Checks on coordinate reference system.
st_crs(province)Coordinate Reference System:
User input: EPSG:32647
wkt:
PROJCRS["WGS 84 / UTM zone 47N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 47N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",99,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand."],
BBOX[0,96,84,102]],
ID["EPSG",32647]]Having a look.
provinceSimple feature collection with 77 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 325178.8 ymin: 620860.6 xmax: 1213656 ymax: 2263241
Projected CRS: WGS 84 / UTM zone 47N
First 10 features:
Shape_Leng Shape_Area ADM1_EN ADM1_TH ADM1_PCODE
1 2.417227 0.13133873 Bangkok กรุงเทพมหานคร TH10
2 1.695100 0.07926199 Samut Prakan สมุทรปราการ TH11
3 1.251111 0.05323766 Nonthaburi นนทบุรี TH12
4 1.884945 0.12698345 Pathum Thani ปทุมธานี TH13
5 3.041716 0.21393797 Phra Nakhon Si Ayutthaya พระนครศรีอยุธยา TH14
6 1.739908 0.07920961 Ang Thong อ่างทอง TH15
7 5.693342 0.54578838 Lop Buri ลพบุรี TH16
8 1.778326 0.06872655 Sing Buri สิงห์บุรี TH17
9 2.896316 0.20907828 Chai Nat ชัยนาท TH18
10 4.766446 0.29208711 Saraburi สระบุรี TH19
ADM1_REF ADM1ALT1EN ADM1ALT2EN ADM1ALT1TH ADM1ALT2TH ADM0_EN ADM0_TH
1 <NA> <NA> <NA> <NA> <NA> Thailand ประเทศไทย
2 <NA> <NA> <NA> <NA> <NA> Thailand ประเทศไทย
3 <NA> <NA> <NA> <NA> <NA> Thailand ประเทศไทย
4 <NA> <NA> <NA> <NA> <NA> Thailand ประเทศไทย
5 <NA> <NA> <NA> <NA> <NA> Thailand ประเทศไทย
6 <NA> <NA> <NA> <NA> <NA> Thailand ประเทศไทย
7 <NA> <NA> <NA> <NA> <NA> Thailand ประเทศไทย
8 <NA> <NA> <NA> <NA> <NA> Thailand ประเทศไทย
9 <NA> <NA> <NA> <NA> <NA> Thailand ประเทศไทย
10 <NA> <NA> <NA> <NA> <NA> Thailand ประเทศไทย
ADM0_PCODE date validOn validTo geometry
1 TH 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((674339.8 15...
2 TH 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((687139.8 15...
3 TH 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((644817.9 15...
4 TH 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((704086 1575...
5 TH 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((662941.6 16...
6 TH 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((643472.8 16...
7 TH 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((751293.3 17...
8 TH 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((647136.1 16...
9 TH 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((620165.4 17...
10 TH 2019-02-18 2022-01-22 -001-11-30 MULTIPOLYGON (((757935.1 16...bkk <- province %>%
filter(ADM1_EN == "Bangkok") # Filter for Bangkok province
print(bkk)Simple feature collection with 1 feature and 16 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 643534.4 ymin: 1492136 xmax: 709544.6 ymax: 1543363
Projected CRS: WGS 84 / UTM zone 47N
Shape_Leng Shape_Area ADM1_EN ADM1_TH ADM1_PCODE ADM1_REF ADM1ALT1EN
1 2.417227 0.1313387 Bangkok กรุงเทพมหานคร TH10 <NA> <NA>
ADM1ALT2EN ADM1ALT1TH ADM1ALT2TH ADM0_EN ADM0_TH ADM0_PCODE date
1 <NA> <NA> <NA> Thailand ประเทศไทย TH 2019-02-18
validOn validTo geometry
1 2022-01-22 -001-11-30 MULTIPOLYGON (((674339.8 15...plot(st_geometry(bkk))
3.0 Geospatial Data wrangling
Check for any duplicated accidents in rdacc_sf. No duplicated accidents events.
any(duplicated(rdacc_sf))[1] FALSEbkk <- st_transform(bkk, crs = st_crs(road)) # Ensure the bkk shapefile has the same CRS
road_bkk <- st_intersection(road, bkk) # Perform the spatial intersection to filter roads in BangkokWarning: attribute variables are assumed to be spatially constant throughout
all geometries#plot(road_bkk) ## don't run me, very longCreate a sf object with only geometry of bkk roads
road_bkk_geom <- st_geometry(road_bkk) #Retain only the geometry of road_bkk
road_bkk_geom_sf <- st_sf(geometry = road_bkk_geom) # Convert the geometry back to a simple feature objectSaving road_bkk_geom_sf
save_path <- "C:/jhui-chen/ISSS626-GAA/Take-Home_Ex/Take-Home_Ex01/data/processed_data/road_bkk_geom_sf.rds" # Specify the file path where you want to save the object
saveRDS(road_bkk_geom_sf, file = save_path) # Save the road_bkk_geom_sf object to the specified fileloading road_bkk_geom_sf back
road_bkk_geom_sf <- readRDS("C:/jhui-chen/ISSS626-GAA/Take-Home_Ex/Take-Home_Ex01/data/processed_data/road_bkk_geom_sf.rds")plot(st_geometry(road_bkk_geom_sf))
4.0 Geospatial data visualisation
4.1 Basic plot of BKK roads and road accidents. This plot provides an overview of the accidents on roads in BKK and serves as a check that the data has been processed correctly. Analysis: There are many roads around BKK area, but it seems that the accidents are clustered around what appears to be expressways and major roads.
plot(st_geometry(road_bkk_geom_sf))
plot(rdacc_sf,add=T,col='red',pch = 19)Warning in plot.sf(rdacc_sf, add = T, col = "red", pch = 19): ignoring all but
the first attribute
4.2 Interactive map of accidents in Bangkok area This plot allows for exploratory data analysis by zooming, panning and inspecting specific locations. Analysis: This helps to show that the accidents are indeed clustered around expressways and major roads.
tmap_mode('view')tmap mode set to interactive viewing tm_shape(rdacc_sf) + # Second shape: accident points
tm_dots(col = "red", size = 0.1) + # Visualize accident points as small red dots
tm_shape(bkk) + # Third shape: Bangkok boundaries
tm_borders(lwd = 2)